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Abstract 

Glaucoma is an ocular disease characterized by progressive retinal ganglion cell (RGC) death caused by axonal injury. 
However, the underlying mechanisms involved in RGC death remain unclear. In this study, we investigated changes in the 
transcriptome profile following axonal injury in mice (C57BL/6) with RNA sequencing (RNA-seq) technology. The experiment 
group underwent an optic nerve crush (ONC) procedure to induce axonal injury in the right eye, and the control group 
underwent a sham procedure. Two days later, we extracted the retinas and performed RNA-seq and a pathway analysis. We 
identified 1 77 differentially expressed genes with RNA-seq, notably the endoplasmic reticulum (ER) stress-related genes Atf3, 
Atf4, Atf5, Chad, Chop, Egrl and Trb3, which were significantly upregulated. The pathway analysis revealed that ATF4 was 
the most significant upstream regulator. The antioxidative response-related genes Hmoxl and Srxnl, as well as the immune 
response-related genes Clqa, Clqb and Clqc, were also significantly upregulated. To our knowledge, this is the first 
reported RNA-seq investigation of the retinal transcriptome and molecular pathways in the early stages after axonal injury. 
Our results indicated that ER stress plays a key role under these conditions. Furthermore, the antioxidative defense and 
immune responses occurred concurrently in the early stages after axonal injury. We believe that our study will lead to a 
better understanding of and insight into the molecular mechanisms underlying RGC death after axonal injury. 
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Introduction 

Glaucoma is a leading cause of blindness worldwide [1]. It is 
characterized by glaucomatous optic neuropathy (GON), and is 
associated with optic nerve degeneration that results in progressive 
visual dysfunction [2] . In glaucoma patients, the number of retinal 
ganglion cells (RGCs) decreases due to axonal degeneration, 
resulting in visual dysfunction. Despite the attempts of many 
clinicians and scientists to identify the molecular mechanisms of 
pathogenesis in glaucoma, they are not yet well understood, 
possibly because of the multifactorial nature of glaucoma [3] . 

High intraocular pressure (lOP) is widely recognized as a major 
risk factor for glaucoma, and treatment to lower lOP is currentiy 
the only method that evidence has shown to prevent the 
progression of the disease [4]. Recently, many varieties of lOP- 
lowering eye drops have become clinically available to treat 
glaucoma. However, substantial reductions in lOP, up to 30%, fail 
to halt the progress of visual dysfunction in some patients, 
particularly those with normal tension glaucoma (NTG) [5]. In 
addition to lOP, risk factors for NTG include age, myopia [6], 
parapapillary atrophy (PPA) [7] and reduced ocular blood flow 
[8] . There is thus a necessity for further investigation of these lOP- 



independent mechanisms, and the development of new neuropro- 
tective drug targets for glaucoma. 

Many recent investigations have led to a growing understanding 
of the underlying mechanism of RGC death in glaucoma, which 
previous studies had found to be induced by axonal injury to the 
lamina cribrosa [3] . However, those studies were mainly designed 
around approaches that focused on only a few pathways [9, 1 0] . In 
order to overcome the heterogeneous and multifactorial nature of 
glaucoma and find novel critical molecular targets for treatment, it 
is necessary to use a global approach (i.e., one including the 
transcriptome and proteome). 

A simple animal model mimicking the pathogenesis of glaucoma 
is a useful tool in investigations of the mechanism of RGC death, 
because standard excisional biopsy is impossible in the case of the 
human retina [11]. In one of the most widely used models, optic 
nerve crush (ONC) is performed in mice to induce axonal injury, 
which is a contributor to the progression of RGC death in 
glaucoma [10,12-14]. Interestingly, in this model the number of 
RGCs is maintained for a short duration after ONC, and 
significant RGC loss is not observed until day 3 [10]. Significant 
axonal damage is known to occur in the retina before visual field 
defects become detectible [15] . It would therefore be very useful to 
develop diagnostic methods and drug targets that functioned in 
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these early stages of glaucoma. Analysis of post-ONC mouse 
retinas in the early stages of axonal injury, before RGC loss (i.e., 
on day 2), may give us valuable insights to help achieve this goal. 

Microarray analysis is a common way to evaluate the expression 
level of large numbers of genes simultaneously. It has also been 
used to evaluate changes after axonal injury both in the retina in 
general and in isolated RGCs [16]. However, microarrays are only 
capable of measuring known transcripts, and do not allow the 
investigation of total genetic changes. By contrast, RNA sequenc- 
ing (RNA-seq) is able to assess complete genes and splice variants, 
with a high degree of reproducibility that matches that of 
microarrays [17]. RNA-seq technology thus has the potential to 
give us very useful, detailed information on the mechanisms of 
disease, as well as unknown pathways and networks of disease, that 
may lead to the discovery of new treatment strategies. 

The purpose of this study was thus to use RNA-serj to 
investigate the molecular mechanisms of damage in the early 
stages of the response to axonal injury, before the onset of RGC 
death. We believe that our study may open new avenues of 
investigation into treatment strategies for axonal damage associ- 
ated with ocular diseases, especially glaucoma. 

Results 

RNA-seq analysis and global gene expression profiles in 
axonal injury 

In order to investigate the transcriptome profile at an early stage 
after axonal injur)-, but before significant RGC loss [10], we 
performed an RNA-seq analysis of mouse retinas harvested 2 days 
after ONC or sham operations. In order to obtain triplicate results, 
three samples were obtained from each group, each sample being 
a combination of material from six unique retinas. Each of these 
samples was sequenced on one lane of the lUumina HiSer-i2000 
platform (lUumina, San Diego, CA). All sequence reads were 
mapped to the reference genome (NCBI37/mm9) with CLC 
Genomics Workbench (version 6.0.1) (CLC Bio, Aarhus, Den- 
mark) [18,19]. The total number of reads per lane was 
approximately 400 million, and the total number of reads per 
sample ranged from 62.9 to 70.3 million. An average 73.8% of 
total reads were mapped in pairs to the reference genome (Table 
SI). Detailed mapping statistics are listed on Table S2. To 
determine the expression level of various genes and compare them 
between samples, we used variable RPKM (reads per kilobase of 
exon per million mapped reads) [20]. To examine the overall 
distribution of gene expression values, we created box plots of 
RPKM expression values with CLC Genomics Workbench 
(Figure lA). Overall RPKM expression values were similar in 
each sample. We excluded genes that did not have a mean RPKM 
> 0.3 in at least one group, in order to remove background noise 
[21]. The number of genes with a mean RPKM > 0.3 in at least 
one group was 13160. These were used for the differential gene 
expression analysis [19]. Fold change (FC) diflferences between the 
mice that underwent ONC and those that underwent a sham 
operation were calculated. The Student's t-test was performed to 
compare the groups with R software (version 3.0.1) [22]. P-values 
were adjusted for multiphcity with the Bioconductor package 
qvalue to control the false discovery rate (FDR) [23]. Differentially 
expressed genes (DEGs) were defined as those with | FC | > 1.5 
and FDR < 0.1 [24]. It is known that the ceUs affected by ONC 
are mainly RGCs. The abundance ratio of RGCs in retinal tissue 
has been reported to very low (less than 0.5%) [25]. Furthermore, 
an approximately 1 .5 to 2-fold increase in gene expression (e.g., in 
Jun, Jmd and Gadd45a) has been reported to be a significant change 
in a previous analysis of changes in the entire retina after optic 



nerve crush [26] . We therefore applied this relatively lower cutoff 
(|FC| > 1.5 and FDR < 0.1) in the current study. We created a 
volcano plot showing DEGs as red dots with the ggplot2 package 
in R software [27] (Figure IB). We also conducted a hierarchical 
clustering analysis of DEGs from all samples with Ward's method 
of Euclidean distances [28], and created a heatmap with the 
heatmap.2 function of the gplots package of R software [29]. The 
results indicated that gene expression was similar in each group 
(Figure IC). 

Differentially expressed genes after ONC 

The triplicate samples from the ONC and sham groups were 
assayed for DEGs, and 177 DEGs (132 up- and 45 downregulated 

genes) were identified (Table S3). The 10 most up- and 
downregulated genes are listed in Table 1. The expression 
changes of known RGC markers and axon regeneration markers 
[16,30,31] are summarized in Table 2. We found that the 
following RGC markers were significantiy downregulated 2 days 
after ONC: Nefh (-2.24-fold), Pou4fl (-1.54-fold), Pou4j2 (-2.40- 
fold), Rhpms (-1.62-fold) and Sncg (-1.77-fold). Interestingly, the 
expression of Thjl and Pou4J3 did not change significantly at this 
time point (day 2). Generally, the expression of Thyl begins to 
decrease 3 days after axotomy in rats [32]. The following axon 
regeneration markers were significantiy upregulated after ONC: 
Gap43 (1.53-fold) and Spnla (23.81-fold). 

A review of the published literature revealed that the following 
sets of endoplasmic reticulum (ER) stress-related genes have been 
shown to be significantly upregulated 2 days after ONC: Atf3, Atf4, 
Atf5, Chad, Ddit3, Egrl, TribS [33-35] (Table 2). 

Significant networks and biological functions after ONC 
revealed by pathway analysis 

To investigate the pathways involved in axoncd injury, the DEG 
dataset was uploaded to Ingenuity Pathway Analysis (IPA, 
Ingenuity Systems, Redwood Cit)', CA) and mapped to the 
Ingenuity Pathways Knowledge Base (IPKB) [36]. The 2 most 
significant networks are shown in Figure 2. Network 1 (Figure 2 A) 
was associated with the "Cell Death and Survival", "Cancer" and 
"Cell Morphology" pathways. Network 2 (Figure 2B) was 
associated with the "Neurological Disease", "Nervous System 
Development and Function" and "Tissue Morphology" pathways. 
Table 3 lists the 5 most significant molecular and cellular 
functions. The most significant biofunction response, according 
to IPA, was for the "Cell Death and Survival" pathway, which 
involved 45 genes (Table S4). 

RT-PCR validation of RNA-seq data 

To validate the RNA-seq findings, we prepared new mouse 
retinas in each group, and performed RT-PCR on these new 
groups of retinas. We selected 14 genes (Spnla, Mmpl2, Soxll, 
Atf3, Tnfrsfl2a, Hmoxl, Plat, Egrl, AtfS, DditS, Jun, Pou4f2, Nejh and 
Pou4fT) involved in the "Cell Death and Survival" pathway, and 
examined changes in their expression with RT-PCR (Table 4). We 
found that results obtained with RT-PCR were similar to those 
obtained with RNA-seq. 

Upstream analysis and global network interactions after 
axonal injury 

In order to investigate molecular network interactions, IPA 
performed an upstream regulator analysis. Table 5 shows the 
transcription factors that IPA predicted to be upstream regulators. 
The most significant was ATF4, but TP53, nuclear factor 
(erythroid-derived 2)-like 2 (NFE2L2) and DNA-damage inducible 
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Figure 1 . Gene expression profiles of the samples. (A) Box plot showing overall RPKM expression values for the ONC and control samples. (B) 
Volcano plot showing differentially expressed genes after axonal injury. For each plot, the X-axis represents log2 FC and the Y-axis represents -log,o 
(P-values). DEGs are shown as red dots. (C) Hierarchical clustering of DEGs after ONC. Red indicates increased expression and green indicates 
decreased expression. DEGs were defined as having absolute FC > 1.5 and a FDR < 0.1. 
doi:10.1371/journal.pone.0093258.g001 
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Table 1. Top 10 upregulated and downregulated genes after ONC. 



Symbol 


Description 


Gene accession 


Fold change 


/'-value 


FDR 


Upregulated 


Sprrla 


Small proline-rich protein 1A 


NM_009264 


23.81 


5.46E-05 


0.026 


MmpU 


IVlatrix metallopeptidase 12 


NM_008605 


17.82 


2.80E-04 


0.045 


Ecell 


Endothelin converting enzyme-like 1 


NM_021306 


15.96 


1 .08E-03 


0.054 


Chad 


ChaC, cation transport regulator-like 1 (E. coli) 


NM_026929 


6.61 


1.44E-05 


0.022 


SoxU 


SRY-box containing gene 11 


NM_009234 


5.92 


2.98E-04 


0.045 


Atf3 


Activating transcription factor 3 


NM_007498 


5.34 


7.10E-05 


0.028 


Lgals3 


Lectin, galactose binding, soluble 3 


NM_001 145953 


4.23 


3.66E-03 


0.071 


Phgdh 


3-phosphoglycerate dehydrogenase 


NM_0 16966 


4.09 


3.36E-04 


0.045 


Cdknla 


Cyclin-dependent kinase inhibitor 1A (P21) 


NM_007669 


4.03 


6.92 E-04 


0.051 


TnfrsfUa 


Tumor necrosis factor receptor superfamily, member 12a 


NM_013749 


3.98 


4.84E-04 


0.048 


Downregulated 


Gm6747 


Predicted gene 6747 


XIV1_003945591 


-3.53 


2.17E-03 


0.061 


Irx2 


Iroquois related homeobox 2 (Drosophila) 


NM_010574 


-2.77 


3.69E-04 


0.045 


Gm7244 


Predicted gene 7244 


NG_019018 


-2.57 


3.35E-03 


0.069 


Rasgrp2 


RAS, guanyl releasing protein 2 


NM_0 11242 


-2.57 


2.13E-04 


0.043 


Tppp3 


Tubulin polymerization-promoting protein family member 3 


NM_026481 


-2.55 


2.27E-05 


0.022 


Kcnd2 


Potassium voltage-gated channel, Shal-related family, memberNM_019697 
2 


-2.47 


6.78E-03 


0.081 


Opn3 


Opsin 3 


NM_0 10098 


-2.44 


3.63E-03 


0.071 


Ctxn3 


Cortexin 3 


NM_001 134697 


-2.42 


3.62E-04 


0.045 


Pou4f2 


POU domain, class 4, transcription factor 2 


NM_1 38944 


-2.40 


7.99E-03 


0.085 


Pvalb 


Parvalbumin 


NM_013645 


-2.37 


5.66E-04 


0.048 


Differences were considered significant when FDR was < 0.1 and |FC| was > 1.5. 
doi:1 0.1 371/journal.pone.C093258.t001 










Table 2. Expression changes in genes associated with RGCs, axon 


regeneration and ER stress after ONC. 






Symbol 


Description 


Gene accession 


Fold change 


/'-value 


FDR 


RGC 


Nefh 


Neurofilament, heavy polypeptide 


NM_0 10904 


-2.24 


3.1 1 E-04 


0.045 


Pou4fl 


POU domain, class 4, transcription factor 1 


Nl\/l_011143 


-1.54 


5.41 E-03 


0.077 


Pou4f2 


POU domain, class 4, transcription factor 2 


NM_1 38944 


-2.40 


7.99E-03 


0.085 


Pou4f3 


POU domain, class 4, transcription factor 3 


NIVl.l 38945 


1.04 


NS 


NS 


Rbpms 


RNA binding protein gene with multiple splicing 


NIV1_019733 


-1.62 


1.43 E-03 


0.056 


Sncg 


Synuclein, gamma 


NIV1_011430 


-1.77 


1 .27E-04 


0.032 


Thyl 


Thymus cell antigen 1, theta 


NIV1_009382 


-1.07 


NS 


NS 


Axon regeneration 


Gap43 


Growth associated protein 43 


NIV1_008083 


1.53 


4.44E-03 


0.073 


Sprria 


Small proline-rich protein 1A 


NIV1_009264 


23.81 


5.46E-05 


0.026 


ER stress 


Atf3 


Activating transcription factor 3 


NIV1_007498 


5.34 


7.10E-05 


0.028 


Atf4 


Activating transcription factor 4 


Nl\/l_009716 


1.61 


5. 65 E-04 


0.048 


AtfS 


Activating transcription factor 5 


NIV1_030693 


2.24 


2.27E-03 


0.062 


Chad 


ChaC, cation transport regulator-like 1 {E. coli) 


NM_026929 


6.61 


1 .44E-05 


0.022 


Ddit3 


DNA-damage inducible transcript 3 


NIV1_007837 


2.15 


1.51 E-03 


0.056 


Egrl 


Early growth response 1 


NM_007913 


2.25 


7.25E-04 


0.051 


Trib3 


Tribbles homolog 3 (Drosophila) 


NIVl.l 75093 


2.89 


2.98E-03 


0.067 


Differences were considered significant when FDR was < 0.1 and |FC| was > 1.5. NS 


= not significant. 









doi:l 0.1 371 /journal.pone.0093258.t002 
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Figure 2. Network analysis of the effect of ONC on gene expression. These post-ONC significance networks were generated by IPA. The 2 
most significant networks are shown. (A) Network 1 was associated with the "Cell Death and Survival", "Cancer" and "Cell Morphology" pathways. (B) 
Network 2 was associated with the "Neurological Disease", "Nervous System Development and Function" and "Tissue IVIorphology" pathways. Red 
indicates upregulated genes, green indicates downregulated genes, and white indicates genes that were not annotated in this RNA-seq result but 
that formed part of the network. 
doi:1 0.1 371 /journal.pone.0093258.g002 



transcript 3 (DDIT3) were also determined to be upstream 
regulators activated after ONC. Data for ATF4, TP53, NFE2L2, 
DDIT3 and the target genes from the dataset were merged to 
create a graphical representation of the network of molecular 
relationships following ONC (Figure 3). 

Discussion 

In this study, we used RNA-seq to examine the global 
transcriptome profile early after axonal injury, before the onset 
of significant RGC death. We identified 177 DEGs including 
previously uninvestigated molecules in ONC. A pathway analysis 
of these DEGs revealed that the most significant biofunction in 
axonal injury was the "Cell Death and Survival" pathway. We 
found that the ATF4-regulated pathway and other sets of ER 
stress-related genes were significantly upregulated, and that 
NFE2L2 was also involved in axonal injury, as an upstream 
regulator. These results point to the critical role that ER stress 
plays in axonal damage-induced RGC death after ONC. 
Furthermore, the molecular mechanism of the response to axonal 
injury also involved antioxidative defense. 

This study relied on an animal model of ocular disease in which 
axonal injury was induced by ONC [10]. Many animal models 
have been used in recent investigations of novel treatment 
strategies for glaucoma, such as neuroprotection. Various methods 
of inducing RGC loss in animals have previously been described, 
including ONC [10,14,37], optic nerve axotomy [32,38,39], 
intravitreal administration of N-methyl-d-aspartate [40^2] or 
Kainic acid [43-45] induction of glutamate excitotoxicity, and 
tumor necrosis factor-a-induced neuroinflammation [46,47]. In 
contrast to other models, ONC and optic nerve axotomy induce 
axonal damage by direct optic nerve injury, which is the main 
pathogenic component leading to RGC death in glaucoma [14]. 
In models using glutamate toxicity, RGC death occurs immedi- 
ately with TUNEL signals detectible within 6 hours after injury 
[48,49]. On the other hand, in models using neuroinflammation, 
RGC death takes a few weeks and only a small number of cells are 
susceptible. Since the number of surviving RGCs did not 
significantiy decrease until 3 days after ONC in mice [10], we 
were able to examine retinas on the second day after ONC and 
investigate the transcriptome profile of axonal injury-induced 
changes before the onset of RGC death. ONC was thus the most 
appropriate model of glaucoma for our study. 



To our knowledge, this is the first report to use RNA-seq 
analysis to investigate the retinal transcriptome profile early after 
axonal injury. Although several researchers have conducted 
microarray analyses of axonal injury [16,50], the molecular 
mechanisms remain unclear. In contrast to RNA-seq, expression 
microarrays have a number of limitations (e.g., reliance on existing 
knowledge about the genome sequence, background noise and 
lower dynamic range). We therefore performed RNA-seq to 
generate a global view of the transcriptome after axonal injury. 

Microarray analysis of rodent RGCs isolated with fluorescence- 
assisted cell sorting (FACS) has already been reported, and 
clarified the mechanism of axon regeneration after optic nerve 
axotomy [16]. Our RNA-seq analysis, by contrast, included cells 
from the entire retina. Since retinal glial cells are also affected by 
axonal damage after ONC [51], we therefore hypothesized that 
RNA-seq analysis of the entire retina would yield information that 
had not been revealed by previous microarray analyses of FACS- 
purified RGCs. 

In axonal injury, RGCs decrease due to retrograde axonal 
degeneration [10]. Several RGC marker genes are known to be 
downregulated in response to axonal injury [30,52]. In the current 
study, Pou4fl (also known as Brn3a) and Pou4J2 (also known as 
Brn3b) were downregulated 2 days after ONC (Table 2). Brn3 is a 
transcriptional factor expressed in the retina [53]. Furthermore, 
Brn3a is known to be a useful RGC marker, which can be used to 
identify and quantify RGCs both in controls and injured retinas 
[30]. In our study, Tliyl, another well-known RGC marker [54], 
did not decrease significantly (Table 2). The loss of Brn3a-positive 
RGCs was detected earlier than the loss of Fluorogold-labeled 
RGCs [30] . Bm3a may therefore be a useful marker for evaluating 
RGC dysfunction in the early stages after ONC. 

Axotomized RGCs are known to show many similar changes in 
gene expression during axon regeneration [55]. We found that 
Sprrla and Gap43, genes that are related to axon regeneration, 
were sigiiificaiidy upregulated in the retina after axonal injury 
(Table 2). These results support previous findings obtained from a 
microarray analysis of FACS-purified RGCs taken from retinas 
subjected to axonal injury [16]. 

ER stress is thought to play an important role in the 
pathogenesis of several neurological disorders [56]. ER stress 
activates three unfolded protein pathways (UPRs) including RNA- 
dependeiit protein kinase (PKR)-like ER kinase (PERK), inositol- 



Table 3. Top 5 nnolecular and cellular functions significantly nnodulated after ONC. 





Category 


P-value 


Number of Molecules 


Cell Death and Survival 


7.45E-07-1 .83E-02 


45 


Cellular Function and Maintenance 


2.81 E-06-1 .83E-02 


41 


Cell-To-Cell Signaling and Interaction 


4E-06-1.83E-02 


40 


Molecular Transport 


5.13E-05-1.17E-02 


42 


Small Molecule Biochemistry 


5.13E-05-1.83E-02 


36 



Significances were calculated with Fisher's exact test. 
Differences were considered significant at the P < 0.05 level. 
doi:1 0.1 371/journal.pone.0093258.t003 
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requiring kinase 1 (IREl) and ATF6. Prolonged ER stress can also 
induce apoptosis [35]. In the current study, the ER stress-related 
genes AtJ3, A0, Atf5, Chad, DditS (also known as G/EBP 
homologous protein (CHOP)), Egrl and Trib3 were significandy 
upregulated 2 days after ONC (Table 2). Furthermore, IPA 
predicted that ATF4 was the most significant upstream regulator 
(Table 5). Under ER stress conditions, ATF4 is induced by 
eukaryotic inactivation factor 2a, downstream of the PERK 
pathway [57]. This suggests that ATF4 is the key upstream 
transcription factor induced by ER stress in the early stages of 
axonal injury. 

IPA also predicted that CHOP was a significant upstream 
regulator (Table 5). CHOP is transactivated by ATF4, leading to 
ER stress-induced apoptosis [58]. Deletion of CHOP has been 
found to promote RGC survival [59]. According to IPA, ATF4 
was an upstream regulator of CHOP (Table 5). This suggests that 
the ATF4-CHOP pathway plays an important role in axonal 
damage-induced RGC death. We also found that Jun was 
significantly upregulated (Table 4). JUN is known to be activated 
by the IREl-JNK pathway under ER stress conditions [60], and 
can induce apoptosis [26]. Furthermore, we found that ER stress 
related-genes such as TrihS and Chad were significantly upregu- 
lated after ONC (Table 2). TRIB3 has been reported to be 
involved in ER stress-induced apoptosis in 293 and Hela cells [35]. 
CHAC 1 is involved in glutathione depletion and ROS generation 
[61] and is a proapoptotic component of the UPR, downstream of 
the ATF4-ATF3-CHOP cascade in primary human aortic 
endothelial cell lines [62] . To our knowledge, the role of TRIB3 
and CHACl has not yet been investigated in the retina. The 
multiple ER stress-related pathways discussed above were 
activated concurrentiy in the retina after ONC. Therefore, a 
network-based approach [63], considering multiple pathways and 
molecules leading to cell death, is likely the best approach to 
treatment aimed at RGC protection after axonal injury, resem- 
bling the approach to photoreceptor protection that targets two 
cell death pathways [64]. 



Additionally, oxidative stress has been implicated in many 
neurodegenerative diseases [65,66] . In the current study, we found 
that the antioxidative response-related genes Hmoxl and Srxnl 
were significantly upregulated 2 days after ONC (Table S3). IPA 
indicated that NFE2L2 was one of the upstream regulators 
activated after ONC, and that the increased expression oi Hmoxl 
and Srxnl was a downstream effect of NFE2L2 activation (Table 5 
and Figure 3). NFE2L2, also known as Nrf2 (NF-E2 related factor 
2), is a potent transcriptional activator and plays a central role in 
inducing the expression of many cytoprotective genes such as 
Hmoxl and Srxnl [67,68]. Its translocation into the nucleus has 
been observed at an early stage after ONC [9]. This study also 
revealed that Cdknla (also known as p21) was significandy 
upregulated, and indicated that it interacts with Nrf2 (Figure 3). 
Cytoplasmic p21 has been reported to enhance axonal regener- 
ation and functional recovery after spinal injury in rats [69]. 
Furthermore, it has been reported that transcriptional activation of 
cytoprotective genes by Nrl2 is potentiated in the presence of p21 
through facilitated stabilization of Nrf2 [70]. In summary, the 
results of our study indicate that the Nrf2-related pathway is 
activated in response to axonal injury, which may be involved in a 
part of the defense mechanism suppressing RGC death and 
promoting axonal regeneration in the early stages of axonal injury. 
Enhancement of the antioxidant response, along with the 
inhibition of ER stress-related pathways (e.g., ATF-CHOP), may 
have a synergistic protective effect against RGC death after axonal 
injury. 

The immune response has been reported to be involved in 
central nervous system (CNS) injury [55]. Our study found that 
Clqa, Clqb and Clqc, components of Clq belonging to the classical 
complement pathway, were significantly upregulated (Table S3). 
Furthermore, Clqa was included in the "Cell Death and Survival" 
pathway according to IPA (Table S4). Clq has been reported to be 
implicated in the pathogenesis of neurodegenerative diseases such 
as Alzheimer's disease [71]. A previous study used a microarray 
analysis to demonstrate that the complement pathway is 
upregulated in the retina 2 days after ONC [72], an observation 



Table 4. RT-PCR validation of the expression of selected genes related to the "Cell Death and Survival" pathway. 





Symbol 


Description 


Gene Accession 


RNA-seq 




RT-PCR 












FC 


/lvalue 


FC 


/'-value 


Sprria 


Small proline-hch protein lA 


NM. 


.009264 


23.81 


5.46E-05 


232.12 


2.79E-03 


MmpU 


Matrix metallopeptidase 12 


NM. 


.008605 


17.82 


2.80E-04 


84.75 


2.27E-04 


SoxU 


SRY-box containing gene 11 


NM. 


.009234 


5.92 


2.98E-04 


10.14 


3.74E-04 


Atf3 


Activating transcription factor 3 


NM. 


.007498 


5.34 


7.10E-05 


5.51 


1.13E-02 


TnfrsfUa 


Tumor necrosis factor receptor superfamily member 
12a 


NM. 


.001161746 


3.98 


4.84E-04 


7.57 


1.27E-03 


Hmoxl 


Heme oxygenase (decycling) 1 


NM. 


.010442 


3.67 


4.82 E-05 


4.50 


1.53E-02 


Plat 


Plasminogen activator, tissue 


NM. 


.008872 


2.26 


5.63E-05 


2.22 


8.82E-03 


Egrl 


Early growth response 1 


NM. 


.007913 


2.25 


7.25E-04 


4.21 


2.15E-02 


AtfS 


Activating transcription factor 5 


NM. 


.030693 


2.24 


2.27E-03 


2.82 


5.50E-05 


DdltS 


DNA-damage inducible transcript 3 


NM. 


.007837 


2.15 


1.51 E-03 


2.37 


8.60E-06 


Jun 


Jun oncogene 


NM. 


.010591 


2.00 


1.14E-04 


2.22 


1 .78E-03 


Pou4f2 


POD domain, class 4, transcription factor 2 


NM. 


.138944 


-2.40 


7.99E-03 


-2.55 


2.34E-06 


Nefh 


Neurofilament, heavy polypeptide 


NM. 


.010904 


-2.24 


3.1 1 E-04 


-2.36 


3.58E-07 


Pou4fl 


POD domain, class 4, transcription factor 1 


NM. 


.011143 


-1.54 


5.41 E-03 


-2.25 


1.55E-03 



Differences between the NC and sham groups were analyzed with the t-test (RT-PCR: n = 6 for each group). Differences were considered significant at the P<0.05 level. 
doi:1 0.1 371/journal.pone.0093258.t004 
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Table 5. Predicted upstream regulators belonging to transcription factors after ONC. 





Name 


Predicted activation 


Activation Z-score 


/'-value of overlap 


Target molecules in dataset 


ATP/1 


Activated 


3.12 


o.yut- 1 O 


AADC aCMC ATP^ /"nU'MIA nniT3 ("ADC 1 ("Al ^/ITUPm 
AAnb, AblMj, A 1 ro, t^UMNIA, UUWD, uAKj, LuALbjS, IVllnrUZ, 

PSATl, SARS, SERPINFl, SHMT2, SLC7A3, SLC7A5, TNFRSF12A, 
TRIB3 


TP53 


Activated 


2.44 


2.83E-05 


ATF3, ATGIO, CIQC, CDKNIA, CLIC4, DUSP1, HMGN2, HMOXl, 
IFI3C, IGFBP3, KRT18, LGALS3, MIV1P9, SERPINB9, SESN2, TMEM43, 
TRIB3 


NFE2L2 


Activated 


2.13 


4.79E-03 


ADCYAPl, ARHGEF3, CELAl, FXYDl, HAX1, HMOXl, PSATl, 
SRXNl 


DDIT3 


Activated 


2.00 


3.20E-06 


ATF3, ITGAM, SARS, TRIB3 



Data were analyzed with Fisher's exact test. Differences were considered significant with a P<0.05 and |Z-score| ^ 2. 
doi:1 0.1 371 /journal.pone.C093258.t005 



that was repeated in our study. The immune system miglit also Conclusion 
play an important role in the pathogenesis of axonal injury. 

We used RNA-seq technology to investigate the entire retinal 
transcriptome profile in the early stages of post-axonal injury. A 
pathway analysis of DEGs indicated that cell death and the 
survival response were induced at an early stage after ONC. ER 




Regulator 
me Receptor 



Figure 3. Interaction networks Involved In axonal injury after ONC. The upstream analysis was performed with IPA. ATF4, TP53, NFE2L2, 
DDIT3 and the target molecules of these upstream regulators were merged for this representation of the interaction networks after ONC. 
doi:1 0.1 371/journal.pone.0093258.g003 
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stress was the main response in axonal injury, inducing many other 
pathways (i.e., RGC marker down regulation, the antioxidative 
response, the immune response, and axon regeneration). Our 
transcriptomic approach to this investigation, which relied on 
RNA-seq, wa.s a powerful and effective method, and allowed us to 
obtain a global view of gene expres.sion changes in tiic retina after 
axonal injury. We believe that our study has provided new insights 
into the molecular mechanisms underlying axonal damage, and 
may help in research aimed at the discovery of new biomarkers 
and therapeutic targets for a variety of ocular diseases. 

Materials and Methods 

Animals 

C57BL/6 mice (male, 12-week-old; SLC, Hamamatsu, Japan) 
were used in this study. The surgical procedures were performed 
under deep anesthesia with intramuscular administration of a 
mixture of ketamine (100 mg/kg) and xylazine (9 mg/kg). All 
animals were maintained and handled in accordance with the 
guidelines of the ARVO Statement for the Use of Animals in 
Ophthalmic and Vision Research and the guidelines from the 
declaration of Helsinki and the Guiding Principles in the Care and 
Use of Animals. All experimental procedures described in the 
present study were approved by the Ethics Committee for Animal 
Experiments at Tohoku University Graduate School of Medicine, 
and were performed according to the National Institutes of Health 
guidelines for the care and use of laboratory animals. 

Induction of axonal injury in mice 

ONC was used to induce axonal injury as pr(;viously described 
[10]. Briefly, the optic nerve was exposed, crushed 2 mm posterior 
to the globe with fine forceps for 10 seconds, and released. A 
fundus examination was used to confirm the appearance of normal 
blood circulation, and antibiotic ointment was applied. The 
operation was similar in the sham group, but after exposure, the 
optic nerve was not crushed. 

RNA extraction 

Two days after surgery, the retinas of the mice were extracted 
and immediately immersed in an RNA stabilization reagent 
(RNase later sample and assay technology; Qiagen, Valencia, CA). 
The retinas were then homogenized in Qiazol (Qiagen) with a 
pesde homogenizer, and total RNA was extracted from the 
homogenized mixture with a miRNeasy mini kit (Qiagen). The 
resulting 48 individual samples (24 from the ONC group and 24 
from the control group) were then assessed with a spectropho- 
tometer to determine their total RNA concentration (NanoDrop 
2000c, Thermo Scientific). 

RNA sequencing 

Thirty-six samples of purified RNA (18 from the ONC and 18 
from the control group) were used for this analysis. In each group, 
fixed quantities of RNA were taken from six samples and 
combined into a single sample, in order to minimize the influence 
of individual variations in the mice. This process yielded three 
combined samples from both the ONC and control groups. The 
quality of these six combined RNA samples was then assessed with 
an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, 
CA). The triplicated ONC and control samples used for the RNA- 
seq analysis had RNA integrity numbers (RIN) ranging from 7.8 to 
8.2 (Table SI). The cDNA library of each sample was prepared 
with Illumina Tru-Seq RNA Sample Prep Kits (lUumina, San 
Diego, CA) for 100 bp paired-end reads, according to the 
manufacturer's instructions. Each of the six cDNA libraries was 



indexed for multiplexing. These six indexed libraries were 
sequenced on one lane of the Illumina Hiseq2000 device. 

Data were recorded in the FASTQ format and then imported to 
CLC Genomics Workbench (version 6.0.1) (CLC Bio, Aarhus, 
Denmark) for analysis [18,19]. All sequence reads were mapped to 
the reference genome (NCBI37/ mm9) with the RNA-seq mapping 
algorithm included in CLC Genomics Workbench. The maximum 
number of mismatches allowed for the mapping was set at 2. To 
estimate gene expression levels, we calculated RPKM with CLC 
Genomics Workbench, as defined by Mortazavi et al. [20], and 
then analyzed differentially expressed g(;ncs (DEGs) in tlu' control 
and ONC samples. AU sequence data have been deposited in the 
Gene Expression Omnibus under the accession number 
GSE55228. 

Quantitative RT-PCR 

Twelve samples of purified RNA (6 from the ONC and 6 from 
the control group) were used for quantitative RT-PCR. Total 
RNA (200 ng per sample) from the samples was first reverse- 
transcribed into cDNA using Superscript III (Invitrogen Life 
Technologies, Carlsbad, CA). Quantitative RT-PCR was then 
performed with a 7500 Fast Real-Time PGR System (Applied 
Biosystems, Foster City, CA) as previously described, with minor 
modifications [73]. For each 20 ^1 reaction the following were 
used: 10 jll TaqMan Fast Universal PGR Master Mix (Applied 
Biosystems, Foster City, CA), 1 [d Taqman probe, 1 |ll template 
DNA, and 8 [il DEPC water. Each sample was run in duplicate in 
each assay. For a relative comparison of gene expression, we 
analyzed the results of the real-time PGR data with the 
comparative Ct method (2 ^^"^t^^ normalized to Gapdh, an 
endogenous control. All Taqman probes used for these reactions 
are listed in Table S5. 

Pathway analysis 

Pathway and global functional analyses were performed with 
IPA software [36,74,75]. The DEG datasets were uploaded to the 
IPA application and mapped to IPKB. Each gene identifier was 
then mapped to its corresponding IPKB. Networks of these genes 
were generated based on their connectivity. The significance of the 
association between the datasets and biofunrtions were measured 
using a ratio of number of genes from the dataset that map to the 
pathway divided by the total number of genes in that pathway. An 
upstream regulator analysis was performed to compare DEGs in 
the datasets to those known to be regulated by a given upstream 
regulator. Based on the concordance between them, an activation 
score was assigned, showing whether a potential transcriptional 
regulator was in an "activated" (z score S 2), "inhibited" (z score 
£ —2), or uncertain state. 

Statistical analysis 

RNA-seq data were analyzed and RPKM was calculated with 
CLC Genomics Workbench [76]. A threshold RPKM value- of 0.3 
has been reported to balance the numbers of false positives and 
false negatives [21,77]. We therefore excluded genes that did not 
have RPKM > 0.3 in at least one group. This yielded 13160 
genes, which we then used in the differential expression analysis. 
P-values were calculated with the Student's t-test and were 
adjusted for multiphcity with the Bioconductor package qvalue 
[78,79]. This software allows for selecting statistically significant 
genes while controlhng the estimated false discovery rate (FDR). 
FDR < 0.1 with I EC I > 1.5 was considered statistically 
significant in the RNA-seq analysis. RT-PCR data were analyzed 
with the Welch's t-test. Statistical analysis of the RNA-seq and 
RT-PCR data was performed with R software (version 3.0.1) [22]. 
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The significance of the pathway analysis was calculated with 
Fisher's exact test in the IPA application. If the Avalues for RT- 
PCR and IPA were less than 0.05, the result was considered 
statistically significant. 

Supporting Information 

Table SI RNA integrity numbers and summary of sequence 

statistics. 
(XLSX) 

Table S2 Detailed mapping statistics. 
PCLSX) 

Table S3 List of DEGs after ONC. 
(XLSX) 
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